05_cba_results

Published

19/04/2023 T16:10

Code
library(here)
library(flextable)
library(plotly)
source(here("02_scripts/03a_fun_calib_td.R"))
source(here("02_scripts/06_interventions_models.R"))

theme_set(theme_minimal())

Scenario 1: Accounts for the effect of Fentanyl contamination starting from 2018 Scenario 2: Scenario 2: Doesn’t account for Fentanyl contamination and has the same probability from Illicit use to OD and from OD to OD death

Code
trace_tbl <- as_tibble(mod_basecase$m_M) %>% 
  mutate(cycle_num = 0:n_cycles) %>% 
  pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
  mutate(mod = "No Interventions") %>% 
  bind_rows(., as_tibble(mod_nalx_1$m_M) %>% 
              mutate(cycle_num = 0:n_cycles) %>% 
              pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
              mutate(mod = "Naloxone")) %>% 
  bind_rows(., as_tibble(mod_ss_1$m_M) %>% 
              mutate(cycle_num = 0:n_cycles) %>% 
              pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
              mutate(mod = "Safer Supply")) %>% 
  bind_rows(., as_tibble(mod_pres_guid_1$m_M) %>% 
              mutate(cycle_num = 0:n_cycles) %>% 
              pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
              mutate(mod = "Prescription Guidelines")) %>% 
  bind_rows(., as_tibble(mod_all_int_1$m_M) %>% 
              mutate(cycle_num = 0:n_cycles) %>% 
              pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
              mutate(mod = "All Interventions"))

trace_tbl$mod <- factor(trace_tbl$mod,
                        levels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"),
                        labels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"))

v_extra_sevbi <- c(mod_basecase$extra_sev_bi, mod_nalx_1$extra_sev_bi,
                   mod_ss_1$extra_sev_bi, mod_pres_guid_1$extra_sev_bi,
                   mod_all_int_1$extra_sev_bi)

v_extra_modbi <- c(mod_basecase$extra_mod_bi, mod_nalx_1$extra_mod_bi,
                   mod_ss_1$extra_mod_bi, mod_pres_guid_1$extra_mod_bi,
                   mod_all_int_1$extra_mod_bi)

v_extra_od <- c(mod_basecase$extra_od_deaths, mod_nalx_1$extra_od_deaths,
                   mod_ss_1$extra_od_deaths, mod_pres_guid_1$extra_od_deaths,
                   mod_all_int_1$extra_od_deaths)

mod_names <- c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions")
mod_names <- factor(mod_names, levels = mod_names, labels = mod_names)

for (j in 1:length(mod_names)){
  
  trace_tbl$count[trace_tbl$state == "BO_OD_DEATH" &
                      trace_tbl$cycle_num == 180 &
                      trace_tbl$mod == mod_names[j]] <- trace_tbl$count[trace_tbl$state == "BO_OD_DEATH" &
                      trace_tbl$cycle_num == 180 &
                      trace_tbl$mod == mod_names[j]] + v_extra_od[j]
  
  trace_tbl$count[trace_tbl$state == "BO_MOD_BI" &
                      trace_tbl$cycle_num == 180  &
                      trace_tbl$mod == mod_names[j]] <- trace_tbl$count[trace_tbl$state == "BO_MOD_BI" &
                      trace_tbl$cycle_num == 180  &
                      trace_tbl$mod == mod_names[j]] + v_extra_modbi[j]
  
  trace_tbl$count[trace_tbl$state == "BO_SEVERE_BI" &
                      trace_tbl$cycle_num == 180  &
                      trace_tbl$mod == mod_names[j]] <- trace_tbl$count[trace_tbl$state == "BO_SEVERE_BI" &
                      trace_tbl$cycle_num == 180  &
                      trace_tbl$mod == mod_names[j]] + v_extra_sevbi[j]
}

trace_tbl$state <- factor(trace_tbl$state,
                              levels = v_state_names,
                              labels = v_state_names)
Code
trace_tbl_2 <- as_tibble(mod_basecase_2$m_M) %>% 
  mutate(cycle_num = 0:n_cycles) %>% 
  pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
  mutate(mod = "No Interventions") %>% 
  bind_rows(., as_tibble(mod_nalx_2$m_M) %>% 
              mutate(cycle_num = 0:n_cycles) %>% 
              pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
              mutate(mod = "Naloxone")) %>% 
  bind_rows(., as_tibble(mod_ss_2$m_M) %>% 
              mutate(cycle_num = 0:n_cycles) %>% 
              pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
              mutate(mod = "Safer Supply")) %>% 
  bind_rows(., as_tibble(mod_pres_guid_2$m_M) %>% 
              mutate(cycle_num = 0:n_cycles) %>% 
              pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
              mutate(mod = "Prescription Guidelines")) %>% 
  bind_rows(., as_tibble(mod_all_int_2$m_M) %>% 
              mutate(cycle_num = 0:n_cycles) %>% 
              pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
              mutate(mod = "All Interventions"))

trace_tbl_2$mod <- factor(trace_tbl_2$mod,
                        levels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"),
                        labels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"))

v_extra_sevbi_2 <- c(mod_basecase_2$extra_sev_bi, mod_nalx_2$extra_sev_bi,
                   mod_ss_2$extra_sev_bi, mod_pres_guid_2$extra_sev_bi,
                   mod_all_int_2$extra_sev_bi)

v_extra_modbi_2 <- c(mod_basecase_2$extra_mod_bi, mod_nalx_2$extra_mod_bi,
                   mod_ss_2$extra_mod_bi, mod_pres_guid_2$extra_mod_bi,
                   mod_all_int_2$extra_mod_bi)

v_extra_od_2 <- c(mod_basecase_2$extra_od_deaths, mod_nalx_2$extra_od_deaths,
                   mod_ss_2$extra_od_deaths, mod_pres_guid_2$extra_od_deaths,
                   mod_all_int_2$extra_od_deaths)

for (j in 1:length(mod_names)){
  
  trace_tbl_2$count[trace_tbl_2$state == "BO_OD_DEATH" &
                      trace_tbl_2$cycle_num == 180 &
                      trace_tbl_2$mod == mod_names[j]] <- trace_tbl_2$count[trace_tbl_2$state == "BO_OD_DEATH" &
                      trace_tbl_2$cycle_num == 180 &
                      trace_tbl_2$mod == mod_names[j]] + v_extra_od_2[j]
  
  trace_tbl_2$count[trace_tbl_2$state == "BO_MOD_BI" &
                      trace_tbl_2$cycle_num == 180  &
                      trace_tbl_2$mod == mod_names[j]] <- trace_tbl_2$count[trace_tbl_2$state == "BO_MOD_BI" &
                      trace_tbl_2$cycle_num == 180  &
                      trace_tbl_2$mod == mod_names[j]] + v_extra_modbi_2[j]

    trace_tbl_2$count[trace_tbl_2$state == "BO_SEVERE_BI" &
                      trace_tbl_2$cycle_num == 180  &
                      trace_tbl_2$mod == mod_names[j]] <- trace_tbl_2$count[trace_tbl_2$state == "BO_SEVERE_BI" &
                      trace_tbl_2$cycle_num == 180  &
                      trace_tbl_2$mod == mod_names[j]] + v_extra_sevbi_2[j]
}

trace_tbl_2$state <- factor(trace_tbl_2$state,
                              levels = v_state_names,
                              labels = v_state_names)

trace_tbl <- trace_tbl %>% 
  mutate(scenario = "Scenario 1") %>% 
  bind_rows(., trace_tbl_2 %>% 
              mutate(scenario = "Scenario 2"))

1 Cumulative OD Deaths

1.0.0.1 Cumulative OD Deaths - Scenario 1
Code
p1 <- ggplot(trace_tbl %>% 
               filter(state == "BO_OD_DEATH") %>% 
               filter(cycle_num != 0) %>% 
               filter((cycle_num + 1) %% 12 == 1) %>% 
               filter(scenario == "Scenario 1") %>% 
               mutate(year = rep(c(2015:2029), length(mod_names))),
             aes(x = year, y = count, colour = mod)) +
               geom_line() +
               geom_point(size = 1) +
               ylab("Cumulative OD Deaths") + xlab("Year") +
               labs(colour = "Intervention")

ggplotly(p1)
1.0.0.2 Cumulative OD Deaths - Scenario 2
Code
p2 <- ggplot(trace_tbl %>% 
               filter(state == "BO_OD_DEATH") %>% 
               filter(cycle_num != 0) %>% 
               filter((cycle_num + 1) %% 12 == 1) %>% 
               filter(scenario == "Scenario 2") %>% 
               mutate(year = rep(c(2015:2029), length(mod_names))),
             aes(x = year, y = count, colour = mod)) +
               geom_line() +
               geom_point(size = 1) +
               ylab("Cumulative OD Deaths") + xlab("Year") +
               labs(colour = "Intervention")

ggplotly(p2)
1.0.0.3 Both scenarios
Code
plotly::ggplotly(ggplot(trace_tbl %>% 
               filter(state == "BO_OD_DEATH") %>% 
               filter(cycle_num != 0) %>% 
               filter((cycle_num + 1) %% 12 == 1) %>% 
               mutate(year = rep(c(2015:2029), length(mod_names)*2)),
             aes(x = year, y = count, colour = mod)) +
               geom_line() +
               geom_point(size = 1) +
               ylab("Cumulative OD Deaths") + xlab("Year") +
               labs(colour = "Intervention") +
               facet_wrap(~scenario) +
               labs(caption = "Scenario 1: Accounts for the effect of Fentanyl contamination starting from 2018 \n
                    Scenario 2: Doesn't account for Fentanyl contamination and has the same probability from \n Illicit use to OD and from OD to OD death"))

2 New OD Deaths per year

Code
inci_od_deaths_basecase <- inci_od_deaths_nalox <- inci_od_deaths_ss <- inci_od_deaths_pg <- inci_od_deaths_all <- inci_od_deaths_basecase_2 <- inci_od_deaths_nalox_2 <- inci_od_deaths_ss_2 <- inci_od_deaths_pg_2 <- inci_od_deaths_all_2 <-  rep(NA, 14)

for (i in 1:14){
  yr <- 2015 + i 
  inci_od_deaths_basecase[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
  inci_od_deaths_basecase_2[i] <- mod_basecase_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_basecase_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
  inci_od_deaths_nalox[i] <- mod_nalx_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_nalx_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
   inci_od_deaths_nalox_2[i] <- mod_nalx_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_nalx_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
  inci_od_deaths_ss[i] <- mod_ss_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_ss_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
   inci_od_deaths_ss_2[i] <- mod_ss_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_ss_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
  inci_od_deaths_pg[i] <- mod_pres_guid_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_pres_guid_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
   inci_od_deaths_pg_2[i] <- mod_pres_guid_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_pres_guid_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
  
  inci_od_deaths_all[i] <- mod_all_int_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_all_int_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
  inci_od_deaths_all_2[i] <- mod_all_int_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_all_int_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
}

inc_od_death_tbl <- data.frame(intervention = c(rep(mod_names, each = 14)),
           year = c(rep(2016:2029, 5)),
           inci_od_deaths = c(inci_od_deaths_basecase,
                              inci_od_deaths_nalox,
                              inci_od_deaths_ss,
                              inci_od_deaths_pg,
                              inci_od_deaths_all),
           scenario = rep("Scenario 1", 14*5)) %>% 
  bind_rows(., data.frame(intervention = c(rep(mod_names, each = 14)),
           year = c(rep(2016:2029, 5)),
           inci_od_deaths = c(inci_od_deaths_basecase_2,
                              inci_od_deaths_nalox_2,
                              inci_od_deaths_ss_2, 
                              inci_od_deaths_pg_2,
                              inci_od_deaths_all_2),
           scenario = rep("Scenario 2", 14*5)))

saveRDS(inc_od_death_tbl, file = "inc_od_death_tbl_mod.RDS")

target_od_deaths <- calib_target_tbl %>% 
  filter(group == "total_od_deaths") %>% 
  select(year, target) %>% 
  rename(inci_od_deaths = target) %>% 
  mutate(intervention = "Target")

saveRDS(target_od_deaths, file = "inc_od_death_tbl_target.RDS")

inc_od_death_tbl$intervention <- factor(inc_od_death_tbl$intervention,
                                        levels = mod_names, labels = mod_names)
2.0.0.1 Scenario 1
Code
p3 <- ggplot(inc_od_death_tbl %>% filter(scenario == "Scenario 1"),
                        aes(x = year, y = inci_od_deaths, color = intervention)) +
               geom_line() +
               geom_point(size = 1) +
                 geom_point(data = target_od_deaths, aes (x = year, y = inci_od_deaths)) +
               ylab("New OD Deaths per year") + xlab("Year") +
               labs(colour = "Intervention") +
                 scale_color_brewer(palette = "Dark2")

ggplotly(p3)
2.0.0.2 Scenario 2
Code
p4 <- ggplot(inc_od_death_tbl %>% filter(scenario == "Scenario 2"),
                        aes(x = year, y = inci_od_deaths, color = intervention)) +
               geom_line() +
               geom_point(size = 1) +
                 geom_point(data = target_od_deaths, aes (x = year, y = inci_od_deaths)) +
               ylab("New OD Deaths per year") + xlab("Year") +
               labs(colour = "Intervention") +
                 scale_color_brewer(palette = "Dark2")

ggplotly(p4)
2.0.0.3 Both scenarios
Code
plotly::ggplotly(ggplot(inc_od_death_tbl,
                        aes(x = year, y = inci_od_deaths, color = intervention)) +
               geom_line() +
               geom_point(size = 1) +
                 geom_point(data = target_od_deaths, aes (x = year, y = inci_od_deaths)) +
               ylab("New OD Deaths per year") + xlab("Year") +
               labs(colour = "Intervention") +
                 scale_color_brewer(palette = "Dark2") +
                 facet_wrap(~scenario))

3 Cumulative Deaths

3.0.0.1 Scenario 1
Code
p5 <- ggplot(trace_tbl %>% 
               filter(state == "BO_DEATH") %>% 
               filter(cycle_num != 0) %>% 
               filter((cycle_num + 1) %% 12 == 1) %>% 
               filter(scenario == "Scenario 1") %>% 
               mutate(year = rep(c(2015:2029), length(mod_names))),
             aes(x = year, y = count, colour = mod)) +
  geom_line() +
  ylab("Cumulative Deaths (excluding OD-deaths)") + xlab("Year") +
    labs(colour = "Intervention") 

ggplotly(p5)
3.0.0.2 Scenario 2
Code
p6 <- ggplot(trace_tbl %>% 
               filter(state == "BO_DEATH") %>% 
               filter(cycle_num != 0) %>% 
               filter((cycle_num + 1) %% 12 == 1) %>% 
               filter(scenario == "Scenario 2") %>% 
               mutate(year = rep(c(2015:2029), length(mod_names))),
             aes(x = year, y = count, colour = mod)) +
  geom_line() +
  ylab("Cumulative Deaths (excluding OD-deaths)") + xlab("Year") +
    labs(colour = "Intervention") 

ggplotly(p6)
3.0.0.3 Both scenarios
Code
plotly::ggplotly(ggplot(trace_tbl %>% 
               filter(state == "BO_DEATH") %>% 
               filter(cycle_num != 0) %>% 
               filter((cycle_num + 1) %% 12 == 1) %>% 
               mutate(year = rep(c(2015:2029), length(mod_names)*2)),
             aes(x = year, y = count, colour = mod)) +
  geom_line() +
  ylab("Cumulative Deaths (excluding OD-deaths)") + xlab("Year") +
    labs(colour = "Intervention") +
    facet_wrap(~scenario))

4 Trace Plots

Code
set.seed(122345)
colors_pal <- Polychrome::createPalette(n_states,
                                        c("#084C61", "#DB504A",
                                                   "#E3B505", "#4F6D7A", 
                                                      "#56A3A6"))
names(colors_pal) <- NULL

plotly::ggplotly(ggplot(trace_tbl %>% 
         filter(grepl("BPO", state)), aes(x = cycle_num, y = count, colour = state)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  ggrepel::geom_text_repel(data = trace_tbl %>% filter(cycle_num == 180) %>% 
         filter(grepl("BPO", state)), 
                  aes(x = cycle_num, y = count, label = state), 
                  check_overlap = T, size = 3) +
  scale_color_manual(values = colors_pal) +
  # facet_wrap(~mod) +
    facet_grid(rows = vars(mod), cols = vars(scenario),
               switch = "y"))
Code
plotly::ggplotly(ggplot(trace_tbl %>% 
         filter(state %in% c("BPO_MISUSE", "BI_ILLICIT", 
                             "BS_OAT_INI", "BS_OAT_MAINT")), aes(x = cycle_num, y = count, colour = state)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  ggrepel::geom_text_repel(data = trace_tbl %>% 
         filter(state %in% c("BPO_MISUSE", "BI_ILLICIT",
                             "BS_OAT_INI", "BS_OAT_MAINT")) %>% filter(cycle_num == 180), 
                  aes(x = cycle_num, y = count, label = state), 
                  check_overlap = T, size = 3) +
  scale_color_manual(values = colors_pal)  +
  # facet_wrap(~mod) +
    facet_grid(rows = vars(mod), cols = vars(scenario),
               switch = "y"))
Code
plotly::ggplotly(ggplot(trace_tbl %>% 
         filter(state %in% c("BS_OAT_SS", "BS_SS")), aes(x = cycle_num, y = count, colour = state)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  ggrepel::geom_text_repel(data = trace_tbl %>% 
         filter(state %in% c("BS_OAT_SS", "BS_SS")) %>% filter(cycle_num == 180), 
                  aes(x = cycle_num, y = count, label = state), 
                  check_overlap = T, size = 3) +
  scale_color_manual(values = colors_pal) +
  # facet_wrap(~mod) +
    facet_grid(rows = vars(mod), cols = vars(scenario),
               switch = "y"))
Code
plotly::ggplotly(ggplot(trace_tbl %>% 
         filter(state %in% c("BO_OD_ILLICIT", "BO_OD_RX",
                             "BO_MOD_BI", "BO_SEVERE_BI")), aes(x = cycle_num, y = count, colour = state)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  ggrepel::geom_text_repel(data = trace_tbl %>% 
         filter(state %in% c("BO_OD_ILLICIT", "BO_OD_RX",
                             "BO_MOD_BI", "BO_SEVERE_BI")) %>% filter(cycle_num == 180), 
                  aes(x = cycle_num, y = count, label = state), 
                  check_overlap = T, size = 3) +
  scale_color_manual(values = colors_pal)  +
  # facet_wrap(~mod) +
    facet_grid(rows = vars(mod), cols = vars(scenario),
               switch = "y"))
Code
plotly::ggplotly(ggplot(trace_tbl %>% 
         filter(state %in% c("BR_ILLICIT", "BR_OAT_INI",
                             "BR_OAT_MAINT", "BR_OAT_SS",
                             "BR_SS", "BR_OD_ILLICIT")), aes(x = cycle_num, y = count, colour = state)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  ggrepel::geom_text_repel(data = trace_tbl %>% 
         filter(state %in% c("BR_ILLICIT", "BR_OAT_INI",
                             "BR_OAT_MAINT", "BR_OAT_SS",
                             "BR_SS", "BR_OD_ILLICIT")) %>% filter(cycle_num == 180), 
                  aes(x = cycle_num, y = count, label = state), 
                  check_overlap = T, size = 3) +
  scale_color_manual(values = colors_pal)  +
  # facet_wrap(~mod) +
    facet_grid(rows = vars(mod), cols = vars(scenario),
               switch = "y"))

5 CBA results

5.0.1 Scenario 1

Code
cost_diff <- c(mod_nalx_1$total_net_present_cost - mod_basecase$total_net_present_cost,
               mod_ss_1$total_net_present_cost - mod_basecase$total_net_present_cost,
               mod_pres_guid_1$total_net_present_cost - mod_basecase$total_net_present_cost,
               mod_all_int_1$total_net_present_cost - mod_basecase$total_net_present_cost)

cost_diff_per <- c(round((cost_diff/c(mod_basecase$total_net_present_cost))*100,2))

deaths_eff_tbl <- trace_tbl %>% 
  filter(scenario == "Scenario 1") %>% 
  filter(cycle_num == 180) %>% 
  select(-cycle_num) %>% 
  pivot_wider(names_from = mod, values_from = count) %>% 
  filter(state %in% c("BO_DEATH", "BO_OD_DEATH")) %>% 
  mutate(nalx_effect = Naloxone - `No Interventions`,
         ss_effect = `Safer Supply` - `No Interventions`,
         pg_effect = `Prescription Guidelines` - `No Interventions`,
         all_effect = `All Interventions` - `No Interventions`,
         nalx_effect_per = round(((Naloxone - `No Interventions`)/`No Interventions`)*100, 2),
         ss_effect_per = round(((`Safer Supply` - `No Interventions`)/`No Interventions`)*100, 2),
         pg_effect_per = round(((`Prescription Guidelines` - `No Interventions`)/`No Interventions`)*100, 2),
         all_effect_per = round(((`All Interventions` - `No Interventions`)/`No Interventions`)*100, 2)) %>% 
  select(-c(Naloxone, `No Interventions`, `Safer Supply`, `Prescription Guidelines`, `All Interventions`))


total_death_oddeaths <- trace_tbl %>% 
  filter(scenario == "Scenario 1") %>%
  filter(cycle_num == 180) %>% 
  select(-cycle_num) %>% 
  pivot_wider(names_from = mod, values_from = count) %>% 
  filter(state %in% c("BO_DEATH", "BO_OD_DEATH"))

od_deaths_diff <- c(deaths_eff_tbl$nalx_effect[1], deaths_eff_tbl$ss_effect[1],
                 deaths_eff_tbl$pg_effect[1], deaths_eff_tbl$all_effect[1])
od_deaths_diff_per <- c(deaths_eff_tbl$nalx_effect_per[1], deaths_eff_tbl$ss_effect_per[1],
                 deaths_eff_tbl$pg_effect_per[1], deaths_eff_tbl$all_effect_per[1])
deaths_diff <- c(deaths_eff_tbl$nalx_effect[2], deaths_eff_tbl$ss_effect[2],
                 deaths_eff_tbl$pg_effect[2], deaths_eff_tbl$all_effect[2])
deaths_diff_per <- c(deaths_eff_tbl$nalx_effect_per[2], deaths_eff_tbl$ss_effect_per[2],
                 deaths_eff_tbl$pg_effect_per[2], deaths_eff_tbl$all_effect_per[2])


saveRDS(cbind.data.frame(cost_diff, cost_diff_per,
                         od_deaths_diff, od_deaths_diff_per,
                         deaths_diff, deaths_diff_per), "point_estiamte.RDS")
saveRDS(total_death_oddeaths, "total_deaths_oddeaths.RDS")

costs <- paste0(round(cost_diff/1000000, 0), " (", cost_diff_per, "%)")
deaths <- paste0(round(deaths_diff, 0), " (", deaths_diff_per, "%)")
oddeaths <- paste0(round(od_deaths_diff, 0), " (", od_deaths_diff_per, "%)")

effects_tbl <- cbind.data.frame(mod_names[-1], costs, deaths, oddeaths)

effects_tbl |> 
  flextable() |> 
  add_header_row(colwidths = c(1, 3),
                      values = c("", "Mean Change Compared to No Intervention")) |>
  set_header_labels(values = list(
                         `mod_names[-1]` = "Interventions",
                         costs = "Discounted Net Present Costs \n in Millions (%)",
                         deaths = "Deaths (%)",
                         oddeaths = "Overdose-related Deaths (%)")) |>
  theme_vanilla() |>
  set_table_properties(layout = "autofit") |>
  align_text_col(align = "center", header = TRUE, footer = TRUE)

Mean Change Compared to No Intervention

Interventions

Discounted Net Present Costs
in Millions (%)

Deaths (%)

Overdose-related Deaths (%)

Naloxone

102 (0.01%)

396 (0.01%)

-6420 (-4.05%)

Safer Supply

4233 (0.53%)

-13956 (-0.31%)

-3138 (-1.98%)

Prescription Guidelines

-26103 (-3.28%)

15252 (0.34%)

-1764 (-1.11%)

All Interventions

-24876 (-3.13%)

14115 (0.31%)

-8543 (-5.39%)

5.0.2 Scenario 2

Code
cost_diff_2 <- c(mod_nalx_2$total_net_present_cost - mod_basecase_2$total_net_present_cost,
               mod_ss_2$total_net_present_cost - mod_basecase_2$total_net_present_cost,
               mod_pres_guid_2$total_net_present_cost - mod_basecase_2$total_net_present_cost,
               mod_all_int_2$total_net_present_cost - mod_basecase_2$total_net_present_cost)

cost_diff_per_2 <- c(round((cost_diff_2/c(mod_basecase_2$total_net_present_cost))*100,2))

deaths_eff_tbl_2 <- trace_tbl %>% 
  filter(scenario == "Scenario 2") %>% 
  filter(cycle_num == 180) %>% 
  select(-cycle_num) %>% 
  pivot_wider(names_from = mod, values_from = count) %>% 
  filter(state %in% c("BO_DEATH", "BO_OD_DEATH")) %>% 
  mutate(nalx_effect = Naloxone - `No Interventions`,
         ss_effect = `Safer Supply` - `No Interventions`,
         pg_effect = `Prescription Guidelines` - `No Interventions`,
         all_effect = `All Interventions` - `No Interventions`,
         nalx_effect_per = round(((Naloxone - `No Interventions`)/`No Interventions`)*100, 2),
         ss_effect_per = round(((`Safer Supply` - `No Interventions`)/`No Interventions`)*100, 2),
         pg_effect_per = round(((`Prescription Guidelines` - `No Interventions`)/`No Interventions`)*100, 2),
         all_effect_per = round(((`All Interventions` - `No Interventions`)/`No Interventions`)*100, 2)) %>% 
  select(-c(Naloxone, `No Interventions`, `Safer Supply`, `Prescription Guidelines`, `All Interventions`))

od_deaths_diff_2 <- c(deaths_eff_tbl_2$nalx_effect[1], deaths_eff_tbl_2$ss_effect[1],
                 deaths_eff_tbl_2$pg_effect[1], deaths_eff_tbl_2$all_effect[1])
od_deaths_diff_per_2 <- c(deaths_eff_tbl_2$nalx_effect_per[1], deaths_eff_tbl_2$ss_effect_per[1],
                 deaths_eff_tbl_2$pg_effect_per[1], deaths_eff_tbl_2$all_effect_per[1])
deaths_diff_2 <- c(deaths_eff_tbl_2$nalx_effect[2], deaths_eff_tbl_2$ss_effect[2],
                 deaths_eff_tbl_2$pg_effect[2], deaths_eff_tbl_2$all_effect[2])
deaths_diff_per_2 <- c(deaths_eff_tbl_2$nalx_effect_per[2], deaths_eff_tbl_2$ss_effect_per[2],
                 deaths_eff_tbl_2$pg_effect_per[2], deaths_eff_tbl_2$all_effect_per[2])

costs_2 <- paste0(round(cost_diff_2/1000000, 0), " (", cost_diff_per_2, "%)")
deaths_2 <- paste0(round(deaths_diff_2, 0), " (", deaths_diff_per_2, "%)")
oddeaths_2 <- paste0(round(od_deaths_diff_2, 0), " (", od_deaths_diff_per_2, "%)")

effects_tbl_2 <- cbind.data.frame(mod_names[-1], costs_2, deaths_2, oddeaths_2)

effects_tbl_2 |> 
  flextable() |> 
  add_header_row(colwidths = c(1, 3),
                      values = c("", "Mean Change Compared to No Intervention")) |>
  set_header_labels(values = list(
                         `mod_names[-1]` = "Interventions",
                         costs_2 = "Discounted Net Present Costs \n in Millions (%)",
                         deaths_2 = "Deaths (%)",
                         oddeaths_2 = "Overdose-related Deaths (%)")) |>
  theme_vanilla() |>
  set_table_properties(layout = "autofit") |>
  align_text_col(align = "center", header = TRUE, footer = TRUE)

Mean Change Compared to No Intervention

Interventions

Discounted Net Present Costs
in Millions (%)

Deaths (%)

Overdose-related Deaths (%)

Naloxone

85 (0.01%)

303 (0.01%)

-4861 (-3.9%)

Safer Supply

4047 (0.51%)

-12743 (-0.28%)

-2435 (-1.95%)

Prescription Guidelines

-26098 (-3.29%)

15308 (0.34%)

-1485 (-1.19%)

All Interventions

-24882 (-3.13%)

14105 (0.31%)

-6607 (-5.3%)

Code
tbl_effect_plot1 <- cbind.data.frame(Intervention = mod_names[-1], cost_diff_per,
                                    deaths_diff_per, od_deaths_diff_per) %>% 
  pivot_longer(2:4, names_to = "group", values_to = "percent_change") %>% 
  mutate(scenario = "Scenario 1") %>% 
  bind_rows(., cbind.data.frame(Intervention = mod_names[-1], cost_diff_per_2,
                                    deaths_diff_per_2, od_deaths_diff_per_2) %>% 
  pivot_longer(2:4, names_to = "group", values_to = "percent_change") %>% 
  mutate(scenario = "Scenario 2")) %>% 
  mutate(group = ifelse(group == "cost_diff_per_2", "cost_diff_per",
                        ifelse(group == "deaths_diff_per_2", "deaths_diff_per",
                               ifelse(group == "od_deaths_diff_per_2", "od_deaths_diff_per", group))))
5.0.2.1 Scenario 1
Code
p7 <- ggplot(data = tbl_effect_plot1 %>% 
               filter(scenario == "Scenario 1"), aes(x = group, y = percent_change))+
  geom_point(aes(colour = Intervention, shape = scenario), size = 1.5) +
  scale_color_brewer(palette = "Dark2") +
  scale_shape_manual(values = c(15, 16)) +
    geom_hline(yintercept = 0) +
  ylab("Percent Change Compared to No Intervention (%)") +
  xlab("") +
  scale_x_discrete(name ="", 
                    labels=c("Costs","Deaths","OD-releated Deaths"))+
  theme_minimal()

ggplotly(p7)
5.0.2.2 Scenario 2
Code
p8 <- ggplot(data = tbl_effect_plot1 %>% 
               filter(scenario == "Scenario 2"), aes(x = group, y = percent_change))+
  geom_point(aes(colour = Intervention, shape = scenario), size = 1.5) +
  scale_color_brewer(palette = "Dark2") +
  scale_shape_manual(values = c(15, 16)) +
    geom_hline(yintercept = 0) +
  ylab("Percent Change Compared to No Intervention (%)") +
  xlab("") +
  scale_x_discrete(name ="", 
                    labels=c("Costs","Deaths","OD-releated Deaths"))+
  theme_minimal()

ggplotly(p8)
5.0.2.3 Both
Code
ggplotly(ggplot(data = tbl_effect_plot1, aes(x = group, y = percent_change))+
  geom_point(aes(colour = Intervention, shape = scenario), size = 1.5) +
  scale_color_brewer(palette = "Dark2") +
  # scale_shape_manual(values = c(15, 16, 17, 18)) +
  scale_shape_manual(values = c(15, 16)) +
    geom_hline(yintercept = 0) +
  ylab("Percent Change Compared to No Intervention (%)") +
  xlab("") +
  scale_x_discrete(name ="", 
                    labels=c("Costs","Deaths","OD-releated Deaths"))+
    theme_bw() +
    facet_wrap(~scenario))
Code
tbl_effect_plot2 <- cbind.data.frame(Intervention = mod_names[-1], cost_diff_per,
                                     od_deaths_diff_per) %>% 
  mutate(scenario = "Scenario 1") %>% 
  bind_rows(., cbind.data.frame(Intervention = mod_names[-1], cost_diff_per_2,
                                     od_deaths_diff_per_2) %>%
              rename(cost_diff_per = "cost_diff_per_2",
                     od_deaths_diff_per = "od_deaths_diff_per_2") %>% 
              mutate(scenario = "Scenario 2"))
5.0.2.4 Scenario 1
Code
p9 <- ggplot(data = tbl_effect_plot2 %>% filter(scenario == "Scenario 1"),
             aes(x = cost_diff_per, y = od_deaths_diff_per))+
  geom_point(aes(color = Intervention, shape = scenario)) +
  scale_color_brewer(palette = "Dark2") +
  scale_shape_manual(values = c(15, 16)) +
    geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  xlab("Mean Change in Costs Compared to No Intervention (%)") +
  ylab("Mean Change in Overdose-related deaths \n Compared to No Intervention (%)") +
  theme_minimal()

ggplotly(p9)
5.0.2.5 Scenario 2
Code
p10 <- ggplot(data = tbl_effect_plot2 %>% filter(scenario == "Scenario 2"),
             aes(x = cost_diff_per, y = od_deaths_diff_per))+
  geom_point(aes(color = Intervention, shape = scenario)) +
  scale_color_brewer(palette = "Dark2") +
  scale_shape_manual(values = c(15, 16)) +
    geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  xlab("Mean Change in Costs Compared to No Intervention (%)") +
  ylab("Mean Change in Overdose-related deaths \n Compared to No Intervention (%)") +
  theme_minimal()

ggplotly(p10)
5.0.2.6 Both
Code
ggplotly(ggplot(data = tbl_effect_plot2, aes(x = cost_diff_per, y = od_deaths_diff_per))+
  geom_point(aes(color = Intervention, shape = scenario)) +
  scale_color_brewer(palette = "Dark2") +
  scale_shape_manual(values = c(15, 16)) +
    geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  xlab("Mean Change in Costs Compared to No Intervention (%)") +
  ylab("Mean Change in Overdose-related deaths \n Compared to No Intervention (%)") +
  theme_bw() +
    facet_wrap(~scenario))

6 Brain injury

7 Prevalence of substance use treatment

8 Prevalence of prescription opioid use

9 Prevalence of illicit use

10 Prevalence of misuse